Crime and Topography in San Francisco

Does topography have an effect on the likelihood of car break-ins? This regression analysis seeks to provide insights into this question using data from San Francisco, which has the ideal geography for this context.

That is one steep street!

Analysis Plan

The specific analysis plan steps to address this research question are found in Table 1 below.

Show code

data.frame("Step" = c(1:9),
           "Description" = c("Identify research question", 
                             "Select key variables (based on existing literature)",
                             "Collect data",
                             "Merge data",
                             "Data description",
                             "Visualize relationships",
                             "Test OLS assumptions",
                             "Conduct regression analysis",
                             "Interpret results")) %>% 
  kable(caption = "Table 1. Analysis Plan") %>%
  kable_paper(full_width = FALSE) %>%
  column_spec(1, bold = T) %>%
  row_spec(0, bold = T, color = "black")

Table 1: Table 1. Analysis Plan

Step	Description
1	Identify research question
2	Select key variables (based on existing literature)
3	Collect data
4	Merge data
5	Data description
6	Visualize relationships
7	Test OLS assumptions
8	Conduct regression analysis
9	Interpret results

1. Research Question

My research question is: "What is the relationship between topography and car break-ins in San Francisco?

Background

This analysis will focus on the link between elevation and “hilliness” in determining the likelihood of car break-ins in San Francisco. Both terrain and motor vehicle crimes are ubiquitous when discussing living or visiting San Francisco. This relationship has been coined by the phrase - “Crime Doesn’t Climb”.

In April 2021, Young-An Kim & James C. Wo published Topography and crime in place: The effects of elevation, slope, and betweenness in San Francisco street segments. This study analyzes the effects of “hilliness” on crime, controlling for socio-demographic characteristics Kim and Wo (2021). However, their analysis did not separate by specific crime categories and instead included all types of crime, including violent, nonviolent, property, etc. My analysis will focus only on car break-ins rather than all crime reports, as I believe that these crimes will have a more significant relationship with topography.

Understanding the relationship between topography and car break-ins can motivate local-level policy decisions, such as implementing measures proven to reduce crime in areas identified as having characteristics more susceptible to break-ins.

2. Select key variables

Table 2 below contains the key measures of interest in our analysis:

Important Variables

Show code

data.frame("Dependent Variable" = c("Crime (specifically car break-ins)", ""),
           "Independent Variables" = c("Elevation", "Slope"), 
           "Control Variable" = c("Median Income","")) %>% 
  kable(col.names = c("Dependent Variable", "Independent Variables", "Control Variable"), caption = "Table 2. Key variables for regression analysis") %>%
  kable_paper(full_width = FALSE) %>%
  column_spec(1, bold = T) %>%
  row_spec(0, bold = T, color = "black")

Table 2: Table 2. Key variables for regression analysis

Dependent Variable	Independent Variables	Control Variable
Crime (specifically car break-ins)	Elevation	Median Income
	Slope

When discussing topography, both elevation and “hilliness,” or slope, are necessary for inclusion. This is because it more accurately captures the effect of local level topography, which is supported by Kim & Wo.

In any econometric analysis, it is vital to control for socio-economic variables. In this case, it could be that higher elevations in the city are more affluent areas, which may have an impact on crime. Thus, we want to include median income as a control variable.

Here is the regression equation:

\[MotorVehicleTheft_i = \beta_0 + \beta_1Elevation_i + \beta_2Slope_i + \beta_3MedianIncome_i + u_i\]

3. Collect data

• Crime
  • Retrieved from the San Francisco Open Data Portal
  • All crime reports from 2018-01-01 to 2021-11-04.
    • Queried to only include “Motor Vehicle Theft”
  • Locations are aggregated to closest street intersections.
• Elevation contours
  • Retrieved from the San Francisco Open Data Portal
  • 5 ft. elevation contours
  • Used directly for elevation and also to derive slope through a geo-processing workflow in QGIS with GRASS (see this video for details on this process).
• Median income
  • Retrieved through the tidycensus package via the US Census Bureau.
  • census Tract detail level from the 5-year 2015-2019 ACS data.

4. Merge Data

The following code chunk demonstrates how the crime, elevation, slope, and income datasets were merged.

Show code

# Find index of nearest contour to each crime
elev <- st_nearest_feature(x = crimes, y = contours)

# Add elevation and binary slope columns
crimes <- crimes %>% 
  st_join(y = census_geom, join = st_within, left = TRUE) %>% 
  mutate(elev = contours[elev,]$elevation) %>% 
  rename(median_income = estimate) %>% 
  select(date_incid, slope, median_income, elev, geometry)

5. Data description

The summary statistics for the data are provided in Table 3 below.

Show code

crimes %>% 
  st_drop_geometry() %>% 
  select(slope, elev, median_income) %>% 
  psych::describe(fast=TRUE) %>% 
  kable(col.names = c("", "Count", "Mean", "SD", "Min", "Max", "Range", "SE"), caption = "Table 3. Summary statistics for Slope, Elevation, and Median Income variables") %>%
  kable_paper(full_width = FALSE) %>%
  column_spec(1, bold = T) %>%
  row_spec(0, bold = T, color = "black") %>%
  row_spec(c(1,2), bold = T, background = "#E5E5E5")

Table 3: Table 3. Summary statistics for Slope, Elevation, and Median Income variables

		Count	Mean	SD	Min	Max	Range	SE
slope	1	24454	4.492966	4.793824	0	68	68	0.0306554
elev	2	24454	129.321379	120.439659	-5	780	785	0.7701841
median_income	3	23764	114039.269567	47501.416992	12340	208425	196085	308.1390882

Show code

slope_box <- ggplot(data = crimes, aes(x = "", y = slope)) +
  geom_boxplot() +
  geom_jitter(aes(color = slope), 
              width = 0.2,
              size=0.4, 
              alpha=0.025, 
              show.legend = FALSE) +
  theme_classic() +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_blank()) +
  labs(x = "Slope (percent)")

elev_box <- ggplot(data = crimes, aes(x = "", y = elev)) +
  geom_boxplot() +
  geom_jitter(aes(color = elev), 
              width = 0.2,
              size=0.4, 
              alpha=0.025, 
              show.legend = FALSE) +
  theme_classic() +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_blank()) +
  labs(x = "Elevation (feet)")

income_box <- ggplot(data = crimes, aes(x = "", y = median_income)) +
  geom_boxplot() +
  geom_jitter(aes(color = median_income), 
              width = 0.2,
              size=0.4, 
              alpha=0.025, 
              show.legend = FALSE) +
  theme_classic() +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_blank()) +
  labs(x = "Median Income (USD)")

elev_box + slope_box + income_box + plot_annotation(title = 'Boxplots of Independent Variables',
                                                    theme = theme(plot.title = element_text(hjust = 0.5)))

6. Visualize relationships

The following plots show the simple relationships between count of car break ins and the three independent variables (elevation, slope, and median income).

Show code

# Group by income
income_summary <- crimes %>% 
  st_drop_geometry() %>% 
  group_by(median_income) %>% 
  summarize(count = n())

income_plot = ggplot(data = income_summary, aes(x = median_income, y = count)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic() +
  labs(title = "Crime and Median Income",
       x = "Median Income (USD)",
       y = "Number of Break-Ins")

# Group by elevation
elev_summary <- crimes %>% 
  st_drop_geometry() %>% 
  group_by(elev) %>% 
  summarize(count = n())

elev_plot <- ggplot(data = elev_summary, aes(x = elev, y = count)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic() +
  labs(title = "Crime and Elevation",
       x = "Elevation (feet)",
       y = "Number of Break-Ins")

# Group by slope
slope_summary <- crimes %>% 
  st_drop_geometry() %>% 
  group_by(slope) %>% 
  summarize(count = n())

slope_plot <- ggplot(data = slope_summary, aes(x = slope, y = count)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic() +
  labs(title = "Crime and Slope",
       x = "Slope (percent)",
       y = "Number of Break-Ins")

elev_plot + (slope_plot / income_plot)

We see there is a negative correlation between elevation and crime. It is possible that this relationship is not linear, so this variable may fit better if we take the natural log. Since slope in this analysis is a binary variable, we simply see that there are more crimes in areas that are not designated high slope. Lastly, we observe a weak negative correlation between median income and car break-ins.

The following plots show crimes over the time period of our data.

Show code

# Group by all three variables
crimes_summary <- crimes %>% 
  st_drop_geometry() %>% 
  group_by(slope, median_income, elev) %>% 
  summarize(count = n())

# Group by date
crimes_ts <- crimes %>% 
  st_drop_geometry() %>% 
  group_by(date_incid) %>% 
  summarize(count = n())

ts_plot <- ggplot(data = crimes_ts, aes(x = date_incid, y = count)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic() +
  labs(title = "Daily Crime 2018-Present",
       x = "Date (daily)",
       y = "Number of Break-Ins")

crimes_monthly <- crimes_ts %>% 
  mutate(month = lubridate::floor_date(date_incid, "month")) %>% 
  group_by(month) %>% 
  summarize(monthly_sum = sum(count)) %>% 
  filter(month != "2021-11-01")

monthly_plot <- ggplot(data = crimes_monthly, aes(x = month, y = monthly_sum)) +
  geom_point() +
  geom_line() + 
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic() +
  labs(title = "Monthly Crime 2018-Present",
       x = "Date (monthly)",
       y = "Number of Break-Ins")

ts_plot + monthly_plot

We can see that there is an upward trend to our crime data over the time period of our data.

Show code

#bbox <- st_bbox(census_geom)

library(tmap)
tmap_mode("view")

tm_shape(crimes)+
  tm_dots("elev")

7. Test OLS Assumptions

Recall our four key assumptions for OLS:

1. The population relationship is linear in parameters with an additive disturbance.

2. Our \(X\) variables are exogenous, i.e., \(\mathop{\boldsymbol{E}}\left[ u \mid X \right] = 0\).

3. The \(X\) variables have variation.

4. The population disturbances \(u_i\) are independently and identically distributed as normal random variables with mean zero \(\left( \mathop{\boldsymbol{E}}\left[ u \right] = 0 \right)\) and variance \(\sigma^2\) (i.e., \(\mathop{\boldsymbol{E}}\left[ u^2 \right] = \sigma^2\))

We assume that assumptions 1 and 2 hold. Assumption 3 holds because we see that the \(x\) variables have variation in the plots above. Below, the residuals are generated from the main regression and are used to assess the three components of assumption #4.

Show code

crime_lm <- crimes_summary %>% 
  add_predictions(model = lm(count ~ elev + slope + median_income, data = crimes_summary)) %>% 
  mutate(res = count - pred)

ggplot(data = crime_lm, aes(x = res)) +
  geom_histogram(binwidth = 5, fill = "darkblue", col = "black") + 
  labs(title = "Residual Plot",
       x = "Residuals",
       y = "Count") +
  theme_classic()

The residuals do not appear to be normally distributed. There is a long right tail with high outliers. This indicates that the regression analysis is systematically under-predicting counts of car break-ins.

Show code

ggplot(data = crime_lm, aes(sample = res)) +
  geom_qq() +
  geom_qq_line() +
  theme_classic() +
  labs(title = "QQ Plot",
       x = "Theoretical Normal Distribution",
       y = "Analysis Distribution")

The qq plot shows that the distribution is mostly normally distributed for values up to approximately 2. This suggests that the model is not normally distributed for extreme high values.

Show code

x <- ggplot(data = crime_lm, aes(x = count, y = res)) +
  geom_point() + theme_classic() +
  labs(title = "Residuals vs Dependent (Break-Ins)",
       x = "Number of Break-Ins",
       y = "Residual Error")

y1 <- ggplot(data = crime_lm, aes(x = elev, y = res)) +
  geom_point() + theme_classic() + 
  labs(title = "Residuals vs. Elevation",
       x = "Elevation (feet)",
       y = "Residual Error")

y2 <- ggplot(data = crime_lm, aes(x = slope, y = res)) +
  geom_point() + theme_classic() + 
  labs(title = "Residuals vs. Slope",
       x = "Slope (percent)",
       y = "Residual Error")

y3 <- ggplot(data = crime_lm, aes(x = median_income, y = res)) +
  geom_point() + theme_classic() + 
  labs(title = "Residuals vs. Income",
       x = "Median Income",
       y = "Residual Error")

(x + y1) / (y2 + y3)

Since assumptions 1, 2, and 3 appear to be satisfied, using OLS in this case is an unbiased estimator. However, residual plots indicate that assumption 4 may not be satisfied and therefore OLS may not be the estimator with lowest variance. Thus, an alternative estimator to OLS may provide estimates of the true population parameters with less variance.

8. Regression Analysis

Model 1

Show code

model <- lm(formula = count ~ elev + slope + median_income, data = crimes_summary)

summ(model, digits = 7) 

Observations	4391 (23 missing obs. deleted)
Dependent variable	count
Type	OLS linear regression

F(3,4387)	151.4730923
R²	0.0938608
Adj. R²	0.0932411

	Est.	S.E.	t val.	p
(Intercept)	10.1472401	0.3589623	28.2682638	0.0000000
elev	-0.0139627	0.0009376	-14.8922804	0.0000000
slope	-0.1468949	0.0213511	-6.8799698	0.0000000
median_income	-0.0000094	0.0000027	-3.4337539	0.0006008
Standard errors: OLS

Model 2: Remove median income variable

Show code

model <- lm(formula = count ~ elev + slope, data = crimes_summary)

summ(model, digits = 7) 

Observations	4414
Dependent variable	count
Type	OLS linear regression

F(2,4411)	213.9377923
R²	0.0884246
Adj. R²	0.0880113

	Est.	S.E.	t val.	p
(Intercept)	9.5252524	0.2311926	41.2005037	0.0000000
elev	-0.0158054	0.0009825	-16.0870466	0.0000000
slope	-0.1572832	0.0232144	-6.7752419	0.0000000
Standard errors: OLS

Model 3: Logarithmic transformation

Show code

model_transformed <- lm(formula = log(count) ~ elev + slope + median_income, data = crimes_summary)

summ(model_transformed, digits = 7)

Observations	4391 (23 missing obs. deleted)
Dependent variable	log(count)
Type	OLS linear regression

F(3,4387)	247.6027651
R²	0.1448023
Adj. R²	0.1442175

	Est.	S.E.	t val.	p
(Intercept)	1.8826251	0.0391027	48.1456352	0.0000000
elev	-0.0019118	0.0001021	-18.7187543	0.0000000
slope	-0.0211530	0.0023258	-9.0948218	0.0000000
median_income	-0.0000014	0.0000003	-4.6866190	0.0000029
Standard errors: OLS

The coefficients above can be interpreted as the effect of the independent variables on car break-ins.

9. Results

• Interpret the simple visualizations

• Interpret the R^2 of the two models

• Interpret the OLS assumptions

• Conclusions: What are our predicted effects of independent variables on dependent

• Conclusions: What are the studies flaws? Biases?

• Conclusions: What are possibilities for improvement

• Add citations for data sources (SF Data Portal and US Census(?))

• Finish in project interpretation

• Review rubric